%{
We compute a dynamic model using the OE algorithm with aggregate shocks

This program codes up a dynamic model of entry and exit of real estate
agents with endogenous building of experience. 

Notes:

* make sure capacity constaint is never reached, otherwise it creates
inconsistencies with some listings lost 
* automated the copmutation of a transition matrix but only for one lag*
* 11/12 added a commission split
* 11/27 add randomness in number of listings
%}

function y=model_passive_learn_dyn_prices_poisson(params, ke, State, kh, ka, T, ks, kw, l_bar, init, psi, alph, bet, S_state, q_vec, q,license_limit,saveas,spec, buyerfix, num_buyers, Vb_guess)
params

commission_split = @(earn)0.1498*earn.^0.1455;
if spec==0
     fprintf(saveas); fprintf('\n');
elseif spec==1
    saveas = [saveas,'_fixed_learn'];
    fprintf(saveas); fprintf('\n');
    disp(['Learning Externality turned off, taken transition matrix from results',saveas,'.mat'])
    eval(['load results',saveas,'.mat P Exit_matrix']); % FIX TRANSITION PROBABILITIES
    P_fixed = P./(1-Exit_matrix);
elseif spec==2
     saveas = [saveas,'_bs'];
    fprintf(saveas); fprintf('\n');
     disp(['Business Stealing Externality turned off, taken transition matrix from results',saveas,'.mat'])
     eval(['load results',saveas,'.mat Ne Ne_add Nb Nb_add EProfit']); % FIX THE NUMBER OF LISTINGS
     Ne_prof = Ne; 
     Ne_add_prof = Ne_add;
     Nb_prof = Nb;
     Nb_add_prof = Nb_add;   
elseif spec==3
     saveas = [saveas,'_flexcommission'];
    fprintf(saveas); fprintf('\n');
     disp(['Flexible Commission framework, no experienced agents charge 0 commission, others charge competitive commission to make sellers indifferent'])
end     
     
%% Paramteres of the model
mu_fc       = params(1);
sigma_fc    = params(2);
c_e         = params(3);
delt        = params(4);

%disc_exp = 0.5;

%% Preliminaries
Same_Th     = repmat(State(2,:)',1,ks)==repmat(State(2,:),ks,1);
E           = repmat(State(1,:)',1,ks);
E_new       = repmat(State(1,:),ks,1);
E_max       = (E==ke)'.*Same_Th;
W           = repmat(State(3,:)',1,ks);
W_new       = repmat(State(3,:),ks,1);
ind         = sub2ind(size(T),max(1,reshape(W,[1,ks^2])),max(1,reshape(W_new,[1,ks^2])));
AgrT        = reshape(T(ind),ks,ks); % aggregate state transitions

%% Solution
%% Step 1: Guess entry rate and a Policy function (and distribution of agents)
M = xlsread('data_for_model/save_calibrationcontall_experience_l1', 'distribution'); 
dist = [];
for i = 1:kh
    matrix = M((i-1)*size(M,1)/3+1:(i-1)*size(M,1)/3+ke,3);
    dist = [dist repmat(matrix,ka,1)'];
end
dist = repmat(dist,1,kw/kh);
r = dist;
lambd = repmat(dist(1),1,kw);
% Start with policy that is closest to data
policy = (ke+1-State(1,:))/(ke+1);
policy = policy';
V_old = icdf('Lognormal',1-policy,mu_fc,sigma_fc);


try
    eval(['load results/results',saveas,'.mat policy dist V_old lambd r;']);
catch
    load latest.mat policy dist V_old lambd r;
end

found1  = 0;
nmain=1;
max_iter = 1000;
gammap = 0.2;
gammal = 0.2;
Np = 2;
Nl = 2.5;
tic;

w_mat = repmat(1:1:kw,ks,1)==(repmat(State(3,:)',1,kw));
num_agents = sum(w_mat.*repmat(dist',1,kw));
% if max(num_agents)>license_limit
%     dist = license_limit*kw/ks.*ones(size(dist));
%     lambd = lambd/3;
% end;

Listings = repmat(State(1,:)-1,ks,1);
temp = factorial(Listings);

%% First step in optimization - get the search parameters right to match: UNCOMMENT when running for the first time. Only needs to be omtimized for the 
% prices and sale probabilities
% ds_optimization(State,psi,kh,ke,ka); % UNCOMMENT when running for the first time
load obj_min_search_one_price_static.mat params_min
[y2,frac_sold_state,sale_price_state,th,et]= ds_urn_matching_one_price_static(params_min,State,psi,kh,ke,ka,0);


%% Second Step
while found1==0 && nmain < max_iter % iterate over policy function
    %% Step 2: Baseline Compute Distribution of agents, sellers, and buyers, under the current policy function 
    found2=0;
    P = zeros(ks); % initialize transition probabilities
    nsecond=1;
    while found2==0 && nsecond < 100

        d_mat = repmat(dist',1,kw);

        %% Distribute listings proportional to experience share
        num_agents = sum(w_mat.*repmat(dist',1,kw));
        total_exp = sum(w_mat.*repmat(dist'.*(State(1,:)-1)',1,kw));
        Ne = delt*S_state(State(4,:))./num_agents(State(3,:)) + (1-delt)*S_state(State(4,:)).*(State(1,:)-1)./total_exp(State(3,:));
        
        %% Compute sale probabilities and profits IMPORTANT DIFFERENCE: 1) When running the baseline calibration, we compute the number of total buyers to match sellers + sale probabilities. 2) When running any counterfactuals we take the total number of buyers as given! And adjust buyer valuation to make sure that that squares with number of sellers and the sale probabilities
        if (strcmp(saveas,'Baseline') && spec == 0) || buyerfix ==0
            num_b = th.*(dist.*Ne); 
            num_buyers = sum(repmat(num_b',1,kw).*w_mat);
            Nb = et.*(delt*num_buyers(State(3,:))./num_agents(State(3,:)) + (1-delt)*num_buyers(State(3,:)).*(State(1,:)-1)./total_exp(State(3,:))); % only count up the successfull transactions

            Buy_p = reshape(et,1,ks)'; 
            Sale_p = reshape(frac_sold_state,1,ks)'; 
            Price_vec = reshape(sale_price_state,1,ks);
 
        else % take the number of buyers as given and re-compute sale probabilities etc. 
            num_sellers_state = dist.*Ne;
            [frac_sold_state,sale_price_state,th,et,Vb] = ds_urn_matching_one_price_static_buyerV(params_min,State,num_buyers,num_sellers_state,kw,w_mat,Vb_guess);
            Vb_guess = Vb;
            num_b = th.*(dist.*Ne);
            Nb = et.*(delt*num_buyers(State(3,:))./num_agents(State(3,:)) + (1-delt)*num_buyers(State(3,:)).*(State(1,:)-1)./total_exp(State(3,:))); % only count up the successfull transactions
            
            Buy_p = reshape(et,1,ks)'; 
            Sale_p = reshape(frac_sold_state,1,ks)'; 
            Price_vec = reshape(sale_price_state,1,ks);
            % Note: Nb are number of buyers served by agents in each
            % experience level/state. It does NOT mean number of buyer in
            % each exp/state market. Buyers sort not based on their agent
            % experience.
        end 
        
        %% compute transition probabilities         
         % exit rate
        Exit_matrix = repmat(policy',ks,1); % exit after transitioning
               
        %% poisson distribution of clients, only last year of experience matters
        lam  = repmat(Ne'+Nb',1,ks); 
        Pr = exp(-lam).*(lam.^Listings)./temp; % to be used for profit calculation
        P = AgrT.*Pr; 
        
        P = P + AgrT.*(E_max.*repmat(1-sum(P,2),1,ks)); % add the remaining probability to overshooting experience  
        if spec ~=1
            P = (1-Exit_matrix).*P;
        elseif spec==1 % no learning externality, transition as before
            P = (1-Exit_matrix).*P_fixed;
        end
        
        % Update the ss distribution following the OE paper
        entry_r = lambd(State(3,:));
        ss_found=0;
        while ss_found == 0
            r_new = r*P+entry_r.*init; %.*(1-policy') if entrants don't exit
            if max(abs((r-r_new)))<10^-3 || max(abs((r-r_new)./r))<10^-3
                ss_found=1;
            end
            r = r_new;
        end
        dist_new = r; %dist_new = r./q_vec;
        dist_new(q_vec==0) = 0;
         
       %  checking        
        if max(abs(dist_new-dist).*(dist>0.01))<1
             found2=1;
        else
            dist= dist + (dist_new-dist)/4;   
            nsecond = nsecond+1;
        end
         

    end    
    
    
    %% Step 3: Given the distribution compute per-period profits and value function in each state

    EProfit =  sum(((Pr/kw).*Listings),2).*((Ne./(Ne+Nb))'.*Sale_p+(Nb./(Ne+Nb))').*psi.*Price_vec'; %% buyer probability of purchase already in Nb; %% buyer probability of purchase already in Nb, /kw because Pr is matrix that replicates experience transition for all of the kw states. Could multiply by ArgT, BUT that one includes the exit probability. Not the cleanest code, but this is correct.
    
    if spec==3
        bet_client = 0.95; 
        % Matrix of lowest sale probability for each of the nine periods
        Mu0 = Sale_p((State(3,:)-1).*ke+1); % pick out the sale probability of the lowest experience agent in the corresponding state

        % what is seller continuation value - working with 0 experience agent with 0 commission to sell the next period etc.        
        EV_S_test = 0.97*Price_vec*0.4/(1-bet_client*(1-0.4)); % first stab at it;        
        EV_S = EV_S_test';
        % do convergeance instead
        converge=0;
        while converge==0
            EV_S_new = Mu0.*Price_vec'.*(1-0.03) + bet_client.*(1-Mu0).*((AgrT/ke)*EV_S);
            max(abs((EV_S_new-EV_S)./EV_S));
            if max(abs((EV_S_new-EV_S)./EV_S))<0.001
                converge=1;
            end
            EV_S = EV_S_new;
        end
                
        % compute commission for each experience level/time period
        Psi_flex = (Sale_p-Mu0).*(Price_vec'*(1-0.03)-bet_client*EV_S)./(Sale_p.*Price_vec');    
        
        EProfit =  sum(((Pr/kw).*Listings),2).*((Ne./(Ne+Nb))'.*Sale_p+(Nb./(Ne+Nb))').*Psi_flex.*Price_vec'; %% buyer probability of purchase already in Nb, /kw because Pr is matrix that replicates experience transition for all of the kw states. Could multiply by ArgT, BUT that one includes the exit probability. Not the cleanest code, but this is correct.  
    end

    % adjust for commission split
    EProfit = commission_split(EProfit).*EProfit;
    
    
    % compute fixed cost
    lintransf = @(V) V; %4.9*V; % why 4.9???SG 2023 changed to 1. I'm so confused.. 
    FC = exp_range('Lognormal',repmat(mu_fc,ks,1),repmat(sigma_fc,ks,1),zeros(ks,1),lintransf(V_old));   
    FC(isnan(FC))=0;
    
    % add profits + continuation value
    V =EProfit + bet*max(0.01,sum(P.*(repmat(V_old',ks,1)-repmat(FC,ks,1)),2));
    policy_new = 1-logncdf(lintransf(V),repmat(mu_fc,ks,1),repmat(sigma_fc,ks,1));
    
    w_mat = repmat(1:1:kw,ks,1)==(repmat(State(3,:)',1,kw));
  
    V_old = V_old + (V - V_old)./((nmain)^gammap+Np);
    policy = 1-logncdf(lintransf(V_old),repmat(mu_fc,ks,1),repmat(sigma_fc,ks,1));
    %% 1: Entrants pay fixed costs as well
    %V_e = (V_old'-FC)*(w_mat.*repmat(init',1,kw))- c_e; % 1 by kw
    %% 2: Entrants do not pay fixed cost, that is in the entry cost
    V_e = (V_old')*(w_mat.*repmat(init',1,kw))- c_e; % 1 by kw
    
    num_agents = sum(repmat(dist',1,kw).*w_mat);
    lambd_new = max(0,lambd.*(V_e+c_e)/c_e); % changed 4 to 1
    lambd_new = lambd_new + 100*(lambd<10^-18).*(V_e>100); %% update from 0 up
    lambd_new = max(0,lambd_new - 2*max(0,num_agents - license_limit));
%     if max(num_agents)>license_limit
%         lambd_new = lambd_new*0.5;
%     end;
    num_agents = sum(repmat(dist',1,kw).*w_mat);
    lambd_new(q==0)=0;
    V_e(q==0)=0;
    % check for convergeance
    if max(abs((policy_new-policy)./policy_new))<0.005 && max(abs((lambd_new-lambd)./lambd_new).*(lambd>1))<0.005 %max(abs((policy_new-policy)./policy_new))<0.001 % SG policy was here twice!!!
        found1=1;
    end
    nmain=nmain+1;
    lambd = max(0,lambd + (lambd_new - lambd)./(nmain^gammal+Nl)); % - max(0,num_agents - license_limit)); % subreact w/e is over the limit
    
     if mod(nmain,49)==0
        [lambd; V_e]
        disp(['Curent iteration: ',num2str(nmain)]);
        disp(['Distribution: ',num2str(dist)]);
        disp(['Curent exit rates: ',num2str(policy')]);
        disp(['Curent entry rates: ',num2str(lambd)]);
        disp(['Curent value of entry: ',num2str(V_e)]);
        disp(['Curent values: ',num2str(V')]);
    end
end
if nmain==max_iter-1
     disp('Maximum iterations');
end

L = S_state;
%stats;

%% This is where the results are saved!
try 
    mkdir results 
end
eval(['save results/results',saveas,';']); %%
%eval(['save results/latest',saveas,'.mat lambd policy dist V_old']);

%% Optimization (not actually done here)
y=1000; % below is the optimization code that is not used. Instead we solve the model on a parameter grid and choose the best calibration in a stata program
%{
%% This evaluates the objective function y - how well the equilibrium with current paramteres match the calibration. To choose the right model we now do this on a grid. So this part of the code is no longer relevant except for the results are saved
M_data = xlsread('../save_calibrationcontall_experience_l1', 'states_exit'); 
exit_rates_data = M_data(:,2:end);
agr_state = M_data(:,1);
w_agr_state_ex = kh*(floor(agr_state./10)-1)+ mod(agr_state,10); % compute the corresponsing w in the model

M_data = xlsread('../save_calibrationcontall_experience_l1', 'states'); 
entry_data = M_data(:,2);
agr_state = M_data(:,1);
w_agr_state_ent = kh*(floor(agr_state./10)-1)+ mod(agr_state,10); % compute the corresponsing w in the model

M_data = xlsread('../save_calibrationcontall_experience_l1', 'distribution_state'); % to get the entry rates
n=size(w_agr_state_ent,1);
temp = size(M_data,1)/n;
M_data = reshape(M_data(:,3),temp,n);
num_newbies = M_data(1,:);
moments_data = [];
moments_model = [];
for t = 1:size(w_agr_state_ex)
    for i = 1:ke
        choose_state = State(1,:)==i  & (State(3,:)==w_agr_state_ex(t));
        exit_rates_model(i) = sum(choose_state.*policy'); % policy_adj if adjusting policy sum or max. this picks out the right one
    end
    moments_data = [moments_data exit_rates_data(t,1:i)];
    moments_model = [moments_model exit_rates_model];
end
entry_model = [];
for t = 1:size(w_agr_state_ent)
      entry_model(t) = lambd(w_agr_state_ent(t))/sum(dist.*(State(3,:)==w_agr_state_ent(t)));
end

moments_data = [moments_data .77 entry_data']; % changes .77 to match A1
b=regress(Ne'+Nb',State(1,:)'); % Nb already only counts the successful transactions
moments_model = [moments_model b entry_model];

weights = ones(size(moments_data));
weights(end-size(entry_data,1)+1-1:end) = 5; % on the entry and on the listing distribution
y = sum(((abs(moments_data - moments_model)./moments_data).*weights).^2);
if max(lambd)<1 % really low entry rate results in division by zero in exit rate calculation and y=nan
    y = c_b;
end

if abs(exit_rates_model(1)-1)<0.0001; y=10000; end
%}

%% Save Results
if spec==0
    save latest.mat lambd policy dist  V_old
elseif spec==1
    save latest_fix_learn.mat lambd policy dist  V_old
elseif spec==2
    save latest_fix_bs.mat lambd policy dist  V_old
elseif spec==3
    save latest_flexcommission.mat lambd policy dist  V_old
end

